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Abstract: Consider a compound Poisson process with jump measure v 
supported by finitely many positive integers. We propose a method for 
estimating v from a single, equidistantly sampled trajectory and develop 
associated statistical procedures. The problem is motivated by the question 
whether nerve cells in the brain exhibit higher-order interactions in their 
firing patterns. According to the neuronal assembly hypothesis (Hebb 
synchronization of action potentials across neurons of different groups is 
considered a signature of assembly activity, but it was found notoriously 
difficult to demonstrate it in recordings of neuronal activity. Our approach 
based on a compound Poisson model allows to detect the presence of joint 
spike events of any order using only population spike count samples, thus 
bypassing both the "curse of dimensionality" and the need to isolate single- 
neuron spike trains in population signals. 

AMS 2000 subject classifications: 62G05, 62E20, 92C20. 
Keywords and phrases: asymptotics, compound Poisson process, empir- 
ical characteristic function, higher-order interactions, jump measure, spike 
train, synchronized activity. 

Received July 2007. 



1. Introduction 

The cell assembly hypothesis (Hebb [13]) postulates dynamically interacting 
groups of neurons as building blocks of representation and processing of infor- 
mation in the neocortex. Synchronization of action potentials across neuronal 
groups was suggested as a potential signature for active assemblies (von der 
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Malsburg In this scenario, the assembly members would exhibit specific 

higher-order correlations. It was shown that single nerve cells are sensitive to 
higher-order correlations in their input, and that they can react to such input in 
a correlation-specific way (Abeles Diesmann et al. Kuhn et al. [11], 17|). 



The statistical procedures for the analysis of higher-order correlations from mul- 



non 



tiple parallel spike train recordings suggested so far (Perkel et al. [23|, Martig 

et al. [iOl, Griin et al. Q, [3, Nakahara and Amari Giitig et al. [uj, 



[I4I) suffer from the "curse of dimensionality": the estimation of the necessary 
parameters (of order 2^ for a population of iV neurons) poses serious problems, 
mainly due to insufficient sample sizes. As a consequence, most attempts to 
directly observe cell assemblies have resorted to pairwise interactions (Brown 
et al. |6(]). However, when relying on second-order approaches one cannot draw 
direct conclusions about the presence of higher-order effects. 

Here we present a novel procedure for the statistical analysis of higher-order 
patterns in massively parallel spike trains that circumvents the above-mentioned 
problem. It involves two elements: (i) We assume a model for neuronal popula- 
tions with a simple descriptive parametrization of higher-order effects. Specifi- 
cally, we employ correlated point processes, where correlations of any order are 
induced by "inserting" appropriate patterns of synchronous spikes, (ii) We base 
our inference on spike counts extracted from a single population (multi-unit) 
spike train, rather than from multiple parallel (single-unit) spike trains. This 
leads to a parsimoniously parametrized univariate estimation problem circum- 
venting the curse of dimensionality, which greatly reduces the demands with 
respect to the size of empirical samples. This strategy also alleviates some of 
the difhculties arising when single-unit spike trains have to be extracted from 
multi-unit recordings; cf. Section [71 

The mathematical framework is as follows. Consider N neurons labelled 
1,2, ... ,N, and let A denote the collection of all non-empty subsets of the 
set = {1, 2, . . . , A^}. Associated with the neurons is a collection of point pro- 
cesses [pp's] that describe their joint activity. For each subset a € A there is 
a pp Xa(t), t > 00 the jump times of which indicate simultaneous firing of 
all neurons i ^ a. Thus at the lowest level of single neurons there are A^ pp's 
Xi{t),t > (i e N), where Xi{t) denotes the number of "single-cell" spikes 
occurring in the i-th neuron during time interval [0, i]. At the level of neuron 
pairs there are (^) pp's X^ijj{t) counting the number of "two-cell" spikes that 
occur simultaneously in neurons i and j up till time t. Next, there are (^) pp's 
A{ij k} (t) counting the number of "triple-cell" spikes occurring simultaneously 
in neurons i,j,k up till time t. And so on. 

The n-cell processes with n > 2 account for possible rt-fold interactions be- 
tween the single neurons. These processes are not assumed to be observable; only 
their superposition is. Our goal is to disentangle this mixture, i.e., to extract 
information about 2-fold, 3-fold, etc. neuron interactions from the aggregate 
activity of the whole population. 

^We represent pp's by their associated (cumulative) counting processes rather than as 
random measures. 
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Concerning the underlying n-cell processes we make the following assump- 
tions. 

(Al) All processes Xa{t), a € A, are statistically independent. 
(A2) For each a ^ A, Xa{t) is a right-continuous Poisson (point) process [ppp], 
with constant jump rate /i(a), and such that Xa{0) — 0. 

The sum activity of the population is described by the process 

Quantitiy Z(t) counts the number of spikes that occur anywhere in the pop- 
ulation up till time t. Note that a (unit) jump of Xa{t) at time t = t, say, 
contributes \a\ spikes (|a| = cardinality of subset a), hence makes Z{t) jump 
upwards by \a\ units at time t = t. (By (Al), (A2), the probability equals zero 
that there is any i > that is a jump time for more than one of the processes 
Xa, a e A.) 

From representation (jl.ip and our assumptions it follows that Z{t), i > 0, is 
a compound Poisson process with jump measure given by a weighted sum of 
Dirac point masses Sn, 

EN . ^ 

i^nSn , where for every n > 1, Vn = / /^(a) (1-2) 

n— 1 ^ — ^a£A, \a\—n 

is the jump rate of the process Yn '■= "^ZaeA \a\=n -^"-^ superposition of all 
n-cell processes. In terms of the latter, one also has the slightly more condensed 
representation 

Z(t) - V nYn{t), t>0. (1.3) 

^ — 

The Poisson process underlying all jumps, Yn, is called the carrier process. 
Its jump rate is J2n — v+. 

The activity of individual neurons can be described by pp's Wi, i € N, where 
Wi represents the spike train at neuron i € N. In terms of the n-cell processes 
Xa these single- unit spike trains are given by Wi = The collection 

{Wi, i G N} of all single- unit pp's provides a multivariate representation of the 
population activity by means of correlated ppp's: each Wi is a Poisson pp with 
intensity Ai = J2a3i correlation (or "interaction") between Wi and 

Wj {i ^ j) results from those n-cell processes Xa (n > 2) which simultaneously 
contribute to Wi and Wj, i.e., correspond to firing patterns a € A including 
both neurons Note, finally, that the superposition of all single-unit processes 
should yield the total activity of the population, which in fact it does: 

= E. E.,. - E. E.,. - E. - ^ (1-4) 

It has to be pointed out that the subsequent developments exclusively refer 
to the process Z and its assumed compound Poisson character. The precise re- 
lations between, and properties of, the single-unit processes Wi and the n-cell 
processes Xa are technically irrelevant as long as they are compatible with the 
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compound Poisson model. They do matter, however, in regard to the interpre- 
tation of the model and the results derived on its basis. The problem to be dealt 
with can now be stated as follows. 



Problem: Suppose one is given the values of process Z{t) at L equidistant time 
points ti = Ih, 1 < I < L, with spacing h > 0. The task is to estimate the jump 
measure z/, or at least, the first M rates Vn [n < M), from these data. 

Suitable estimates of the Vn could be used to answer specific questions in re- 
gard to neural networks, such as: Are pairwise correlations sufficient to describe 
the network activity? Do there exist correlations of order > 10, say? Etc. 

Note that continuous observation of process Z{t) would allow registration 
of every single jump, hence make the problem quite simple. In practice, one 
depends on discrete sampling, in which case it is unknown how many jumps 
have contributed to the increment across an interval (unless the increment is 
< 1). The information left after discretization is contained in the increments 
of Z{t) across the bins (i;_i,ii]. These bin counts represent poissonized ver- 
sions of the jump measure, which is why we call our goal de-Poissonization. 
De-Poissonization as a powerful analytic method has found numerous applica- 
tions in, e.g., information theory, the analysis of algorithms, and combinatorics 
(Jacquet and Szpankowski Borodin et al. [l]), but we are not aware of ap- 
plications similar to those presented here. For a techni cally related, continuous 
variant of our estimation problem see Jongbloed et al. [l5|. 

The paper is organized as follows. In Section[2]we introduce the rate estimates 
using Fourier inversion. Their basic asymptotic properties are derived in Section 
[3l A more detailed study of the conditions for the validity of the underlying ap- 
proximations is deferred to Appendix [A] The conncection between our approach 
and empirical de-Poissonization is explained in Section Asymptotics-based 
statistical procedures are outlined in Section [5] and illustrated by Monte Carlo 
simulations in Section [6l The paper concludes with a discussion of some neu- 
robiological and statistical aspects of our approach, addressing particularly the 
possible consequences of violations of the compound Poisson assumption and 
possible further developments. 



2. Estimating the rates Un — a Fourier-based proposal 



Given the bin width > 0, let Zi, . . . , denote the increments of process Z{t) 
across the successive sampled time points, Zi = Z{lh) — Z{{1~ l)h). (To ease no- 
tation, we suppress the parameter h wherever feasible.) The Zi are independent 
and identically distributed random variables with the common characteristic 
fimction 



jh{0) = exp 



h 



AxO 



1) v{dx) 



exp 



/i^"^ ^^„(e»^-l) 

^ — ^n— 1 



(2.1) 



Our (first) approach to the estimation of the z^„ depends on the simple ob- 
servation that the exponent of (|2.ip is essentially a Fourier series from which 
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its coefficients, the rates i>m can be recovered by Fourier inversion. With 

J2n=l ^^^^ 



1 

2^ 



V z/„e™^ = i.+ + ;i-ilog7„(0), (2.2) 

^ — ^n— 1 

so by inverting (|2.2p we get back the coefficients of interest, 

{u+ + h-Hog-fh{O))e-"'''d0 (2.3) 

= ^ r h-'\ogji,{e)c-'"'de (n = l,...,M). 

Eq. (|2.3p suggests to estimate the f„ by substituting the unkown character- 
istic function 7^ by the empirical characteristic function [ecf] of the observed 
increments, 

so that the Vn are estimated as fohows, 

^n^^l h-Hog%{9)e-"''d9 (n = l,...,M). (2.4) 

Note that the number of neurons, TV, does not enter p.4|) . hence need not be 
known. The exponential form of 7/1 immediately yields the expression for the 
log characteristic function. Matters are more involved with the empirical ver- 
sion 7/1 whose (continuously defined) complex logarithm may traverse different 
branches. Although such happens only with small probability if L is large (by 
the consistency of the ecf as L f 00, uniformly in 6* £ [— tt, tt]; see e.g. Prakasa 



Rao [2J, Ch. 8.3]), related difficulties may well occur in finite samples. Ways 
of fixing them will develop from an alternative view of our estimation problem 
outlined in Section 3) It has to be emphasized that the estimates can assume 
negative values. We do not attempt to avoid this drawback, and rather consider 
such a case as indicative of a small value of the corresponding rate parameter. 

Linear functionals A = ^„ Cn Vn of the Vm with finitely many real constants 
Cn, are straightforwardly estimated by substituting Un by their estimates, 



/ h-Hog^h{e)c{e)de, (2.5) 

271- 



where in the last expression C{9) — Cfc e ^'"'^ denotes the Fourier series of 
the Cfc (which appears when interchanging summation and integration) . Of par- 
ticular interest in our context are functionals of the form = Tlm=m ^n- Here 
Am represents the total activity of all n-cell processes of orders n > m, and the 
question is whether Am = for all m greater than some mo (and if so, what is 
mo). 

Since the number N of neurons generally is unknown, it would be desir- 
able to replace Am by the infinite tail sum pm = SJ^m'^"' with the under- 
standing that Vn — for all but finitely many n. This is easily achieved on 
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noting that the Fourier series associated with Am is J2n=m e~'"* = j-g-imff _ 
g-i(jv+i)ei-j^^2 — c^*^) and using the Riemann-Lebesgue lemma, according to 
which hmAT^oo /_ f{0) e~*^^ d6 — for every integrable function / S [— tt, tt]. 
Thus 

Pm = lim — / h-'^\og-fh{9) —r: d0 



N^oo 2n 1^ 1 — e 



1 Q-im9 

= — h-Hogjh{e)R,n{e)d0 where = ^; (2.6) 

and at least on any event where (log 7/1 (61))/ (1 — e~'^) G [— tt, tt] it makes 
sense to estimate p,„ by the expression 



Prn = 2^ / log % (9) R,n (9) d9 . 



(2.7) 



The integrability condition holds in fact with probability tending to 1 as L t oo- 
For the moment it suffices to note that the pole of Rm at 6* = is counter- 
balanced by the zero at = of the log (empirical) characteristic function. 

3. Asymptotics 

The basic properties of the proposed estimates that are required for, e.g., test 
construction can be derived applying the usual methods of asymptotic statistics 
(e.g., Prakasa Rao [24]). This is done here in a heuristic fashion at first. Precise 
conditions under which the approximations can be expected to be valid are 
provided in Appendix |^ 

By the law of large numbers, the ecf 7/1(6') is close to its expected value, 

^h{9), with high probability as L becomes large. Let £,h{9) — ^ 1- The 

Taylor expansion 

h'Hog%{9) - h-^\og^h{e) + h-Hog{i + ue)) 

= Ell ^ ^) + (^''^^^ - \^^^^^'' + ■ • ■) (3-1) 

inserted into (|2.4p . with only the first term kept, gives the stochastic approxi- 
mation 

1 r-^ 

(3.2) 



From (j3.2[) one readily obtains the approximate (normal) distribution of the esti- 
mates Vn on calculating moments of [/„ and applying the central limit theorem. 
Clearly is Hermitean, i.e., £,h{9) = ^h{—9) for real 9 (the bar denotes complex 
conjugation), and E^h{9) = 0. Therefore EUn — 0, and by the independence of 



W. Ehm et al./ Empirical de- Poissonization 



479 



the random variables Zi one has for real 6i, 9- 



2 



1 / Qi{ei-e2)Zi \ 



[hL) 

1 / lhi.ei-92) 



= T-^Th{ex,~e2), 



1 



where T = hL denotes "observation length," and Vh the kernel 

h V 7/1 (^1)7/1(^2) / 

Hence, using (|3.2p along with the variable substition —O2 ^ O2 in the integration, 
we find that the asymptotic covariance of estimates Dm , Dn is given by 



= T ^ f^m.n ■ (3-5) 

Slightly more generally, the asymptotic covariance of the estimates Xb, Ac 
of two linear functionals corresponding to finite vectors of constants (bk) = 
b, (cfc) = c is given by the bilinear form 

T-i 6'l]c = rifej &fe ci. 

Alternatively, one may pass to the Fourier series 

B{9) = 6fce"*'=^ C{9) = c^e-^'^' 

and calculate the asymptotic covariance as 

ascov(Ab, Ac) = -|^y" J Vh{9i,92)B{9^)C{92)d9id92. (3.6) 

The latter procedure also applies to certain infinite compounds such as the tail 
sums pm, in which case the Fourier series B, C become i?,„, i?„, respectively; 
see (|2.7p . (Note that F/j (01,6*2) vanishes along the coordinate axes 9i = and 
02 = 0, cancelling the poles of Rm{9i) and i?„(02)-) 

Some insight into the structure of the (co-)variances can be gained by letting 
the bin width h tend to zero. Then by and (PTjl 

N 
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When inserted into p.Sp . this gives the approximation 

^m,n = Vn Sm^n + 0{h) 

for (T times) the asymptotic covariances of the j7„ , which are thus uncorrelated 
to first order, with variance Vn/T . These approximative variances are optimal 
in the sense that they are identical with those applying to the case where each 
n-cell ppp is directly observable. This is in fact little surprising considering that 
if hv^ is small, then any bin contains at most one jump of the carrier process 
with high probability. Note, finally, that r2n,„ = Un vanishes if Vn = 0, indicating 
that Vn tends to zero at a faster rate than Op(r^^/^) in this case. Of course, 
these conclusions apply strictly only in the limit h [Q. 

Let us end this section by indicating the kind of conditions that turn out, in 
Appendix [XI to be decisive for the validity of the proposed normal approxima- 
tions: T should he "large, " h and /T "small, " and hv^ is a critical quantity 
that should be "moderate" at most. 



4. Connection v^rith analytic de-Poissonization 

The estimates (j2.4p can also be derived in a different manner using analytic 
de-Poissonization. This approach provides complementary information about 
our estimates that affords alternative ways of calculating the relevant quantities 
(other than by numerical integration), or even different estimation methods. 

To achieve the necessary generality we will temporarily modify, and slightly 
abuse, the notation used so far. Let us first assume that 

(t>{w) = Ew^ = y^°° pkw'' 

is the generating function of some, for the moment arbitrary, discrete random 
variable Z > with P [Z = k] — pk {k > 0). Suppose that (j) is analytic in 
a neighborhood W of the closed unit disc in the complex plane, and that it 
has no zeros within W. Then log0, too, is analytic in W, and an application of 
Cauchy's integral formulgfl gives (on picking the main branch, where log 1 = 0) 

1 f log^M,^^i,g^(o)^iog,,. (4.1) 



27Ti J\w\ = l 

More generally we have for any n > 

(log0)(")(O) 1 f log (l^{w) ^ 1 



n 



! 27ri 



r i^^,_iri,g^(,.) (4.2) 

Jm=i 27r7_^ 



Let Pn denote the expression on the right-hand side of (|4.2p . Since 0("^(O)/n! 
Pn, (|4.2p then yields the relations 



r, 1 n Pi n P2 ^ f Pl\ n P3 P2 Pi , ^ f Pi 
Po = logpo, /?! = — , /32 = o ~ ' = — 

Po Po 2 \poJ po Po 3 \po 



(4.3) 



For Cauchy's formula and other related facts see Ahlfors 0. 
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Returning to our setting we now let Z stand for the generic specimen of the 
increments Zi. Then 0(e*^) — 7/i(0), and a comparison with (|2.ip to (|2.3p shows 
that the right-hand side of (|4.2p becomes hi/n for n > 1, and —hv+ for n = 0. 
Together with the evaluation (|4.3[) this implies that i/ie raies are expressible 
as functions of the probabilities pk , < k < n, 

hiy+ = -\ogpo, hvi^—, hv2 = —~\i—] , .... (4.4) 

Po Po 2 \poJ 

Moreover, since v+ = pi we get for the tail sums pm 

Pi 

hpi = - logpo , hp2 hpi - hui = -\ogpo , .... (4.5) 

Po 

Taylor expansion of log4>{w) about w = is another way of arriving at these 
relations (Dayley and Vere-Jones 0, Sect. 5.2]); however, the use of Cauchy's 
formula fits better with the Fourier approach. The relations converse to (I4.4p . 

Po = e-^''+, pi=e-''-^hvi, p^^e-'^^^ {{hvif/2 + hv2). ■■■ (4.6) 

supply a basis for alternative estimation approaches to be explored elsewhere. 

A key point of the preceding considerations is that the relations (|4.ip to (|4.5p 
do not only apply to the "true" quantities Vn, pk, but likewise to the estimates 
9„, the empirical frequencies pk = L^^ \{l : Zi = A;}|, and the empirical charac- 
teristic function 'yh{t) — ^k>oPk 6***' — provided the latter has winding number 
no = w.r.t. the origin0 For in this case 7h(t) = 7/i(e**) can be continued to an 
analytic function 7h(w) — X^fePfe polynomial, in fact, since with probabil- 

ity 1 only finitely many pk are nonzero) that is zero-free within the unit disc, 
the latter by the argument principle. Thus, if the winding number is zero, the 
arguments that led to Eq. (14. 2p are valid for the ecf, too, whence 

^^^^i^^ = log7.(e-) e--" dO = (n > 0). (4.7) 

The important conclusion (from the second equality) is that the estimates 9„ 
defined by {2.4-^ and |/^. 7| ), respectively, are identical if uq = 0. 
By evaluating the left-hand side of (|4.7p one obtains as above 

hd+ = - logpo , hdi = pi/po , hV2 = P2/P0 - ^ {Pi/PoY , ■■■ (4-8) 

if fiQ = 0, showing that the rate estimates can be directly calculated from the 
histogram of the bin counts in this case. Explicit expressions for asymptotic 
(co-)variances can be derived in a similar way; cf. Appendix |B1 Unfortunately, 
formulae quickly become complicated, and a computer algebra system will be 
required to handle all but the most simple cases. 

^The curve [— tt, tt] 9 1— > (9) traces out a closed loop within the unit disc of the complex 
plane. Loosely speaking, the winding number w.r.t. counts how often that loop circles around 
(counter-clockwise). For the precise definition see any text book about complex analysis, 
e.g. Ahlfors 0. 
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5. Tests 

Tests of null-hypotheses pertaming to the parameters i^k can be constructed on 
the basis of our estimates and (co-)variance approximations under the assump- 
tion that the estimates are consistent and approximately normally distributed. 
Practical implementation of any such test requires an estimate of the asymptotic 
(co-)variances, that is, of the kernel Tfi{Oi, 62)] the other quantities are known. 
One obvious possibility is to replace the "true" characteristic function 7/1 by its 
empirical version 7/1 everywhere in (|3.4|) . For an alternative that is expected to 
be more stable statistically, note that by (|3.4p and (|2.ip the kernel Th{0\,92) 
admits the explicit representation 



llexp 



n=l 



- 1 , (5.1) 



which allows us to estimate Th{0i,92) by substituting estimates Vn for the un- 
knowns Vn- The P„ can assume negative values, with potentially unpleasant 
consequences. To avoid such, we plug in positive parts 9+ instead (x+ = a; if 
x > 0, and x+ = otherwise). In practice the sum in (|5.ip has to be truncated, 
and the proposed estimate of Vh{Oi, O2) becomes 



%) = \\ exp 



- 1 (5.2) 



with some suitably chosen number K. For example, if a null-hypothesis implies 
Vk = Q for k > m, one may take K = m. This choice is not mandatory, however. 
One also could set = Af if i^i , . . . , vm are the rates that are being estimated^ 
Whatever the concrete choice of F/j, let T~^r2„_„, T~^Yim,m ■■■ denote the 
estimated covariances obtained by replacing F^ in l|3.5p . (jB.ip . ... by F/j. 

Standard ANOVA type tests impose linear restrictions on the rates such as 
Tia: "Av — 0," where v denotes the column vector with entries vi, . . . , I'm and 
A is some qxM matrix of full rank q. The corresponding Wald type test statistic 
W — T {AV)' {AQA')^^ A9 is approximately x^-distributed with q degrees of 
freedom under Ti.Q. Related two-sample analogs may be used to test for equality 
of V in different subjects, or under different conditions in the same subject. 

Other tests relevant to our goals pertain to the tail sums pm- Let F denote the 
maximal index k such that Vk > 0. Clearly then, pm > for every m < V and 
Pm = for every to > 17. This suggests to consider null-hypotheses of the form 
T^o- ^^Pm — y nil < 'Ti < TO2," for some fixed pair toi, TO2 such that 2 < mi < 
TO2. An associated test would reject Ho if V ~ maxmi<m<m2 exceeds some 
critical value, where Vm = {T/Tim,mY^^ Pm- An approximate null-distribution 
for V is difficult to derive theoretically (except for the trivial case toi ~ TO2), 
but may be estimable by means of the bootstrap method. Alternatively, one 
may test each single hypothesis "p™ = 0" observing that Vm is approximately 



* Ultimately, the choice of K should depend on considerations related to the thereby 
achieved validity and efficiency of the resulting test. 
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standard normally distributed (if p,„ = 0), and finally apply some test level 
correction. Most useful in applications may be a simple screening procedure: 
plotting the standardized statistics Vm (or the corresponding p-values) against 
m could suffice to reveal the major features of interest. 



6. Simulation results 



In the following we present simulation results for several parameter configura- 
tions, including one where the Fourier inversion method fails to perform reliably. 
For Example 1 we chose {i^i, . . . , v^) = (40, 10, 4, 3, 1) and i^k = for all k > 5. 
Observation length was T = 30, and bin width h — 0.02. (All times are specified 
in seconds, all rates are given in events/second.) For Example 2 we put i^i = 150 
and i/7 = 7, and Vn = Q otherwise. In this case we had T = 60 and h = 0.005. For 
each example we generated 50 realizations of the respective compound Poisson 
process by Monte Carlo simulation. 



Example 1 
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Example 2 
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Fig 1. Simulated examples. 

Graphical presentation of a numerical simulation of Example 1 (a, c) and Example 2 
{b,d), see text for parameters. a,b, Raster plots, horizontal ticks mark the occurrence 
of a spike in the indicated neuron, some spikes are aligned into spike patterns (top). 
Discrete samples of the compound counting process (bin counts, bottom) shown for 
2 seconds simulation time. c,d, Histogram of bin counts for a full-length simulation. 



In order to illustrate the connection with the usual multivariate approach we 
split each process into a multivariate Poisson process representing a hypothetical 
population of neurons (cf. the paragraph preceding Eq. (|1.4I) ). The populations 
comprised = 30 and N — 20 neurons, respectively, in the two examples 
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shown, yielding average single-neuron firing rates of about 3 Hz and 10 Hz. For 
simplicity, we chose a maximum entropy representation where all neurons have 
the same likelihood to participate in a particular pattern. That is, we generated 
independent ppp's Yn with intensities Vn and assigned the n spikes occurring at 
each jump of Yn randomly to one of the (^) n-element subsets a of the neurons. 
Recall, however, that the single spike trains merely are shown for illustration: 
our statistical analyses only depend on the bin counts of the compound process. 
These are displayed below the single spike trains in Figure [TJ 

For each simulated realization we estimated P„ and p„ for n = 1,...,12 
using the Fourier method, with associated covariance matrix estimates obtained 
according to (|5.2p . As an illustration of hypothesis testing we calculated test 
power profiles /?„ as follows: for each n = 1, . . . , 12 we determined /?„ as the 
relative frequency among the 50 simulations in which the test statistic Vn = 
(r/I]„^„)-^/^ Pn exceeded the value 2, roughly corresponding to a one-sided 5% 
level test of the null-hypothesis "p„ = 0" . The estimated rate profiles reproduce 
the "true" profiles fairly well (Figure [21 left-hand column) . The sharp decrease 
of the power profile in Figure [2ji) properly reflects the vanishing of all t'„ for 
n > 7. Such a clear separation cannot be expected in Example 1 where the 
transition from the non-zero to the zero rates is rather smooth; see Figure [2j:). 



Example 1 Example 2 



c 




Fig 2. De- Poissonization of sample data. 
Fourier-based de- Poissonization applied to the examples described in the text. a,b, 
Shown are the true {u„, black) and estimated [un, blue) rates, plotted against jump 
height n ("Amplitude"), for each of the 50 realizations. c,d, Test power profiles /3„ 
obtained by Monte Carlo simulation, corresponding to null-hypotheses "p„ = 0" . 



Simulations sometimes exhibit a peculiar phenomenon: the 9„ tend to form 
clusters. While in benign cases there is only one cluster concentrated at the 



W. Ehm et al./ Empirical de- Poissonization 



485 



true rates, increase of induces splits of the estimated rates into two or 
more distinct clusters located at totally wrong rate profiles. The reason for 
this phenomenon is that the ecf has a nonzero winding number no in such a 
case, making the log-ecf end up in a different branch. This gives us a simple 
diagnostics: the phase angle of 7/1 (0), log7?i(6'), then changes by a nonzero 
amount 27rno along the loop [— 7r,7r] 3 9 i—f "fh{0). In particular, 31og7ft(7r) — 
3 log 7/1(0) equals a nonzero multiple of tt in such a case (and is zero otherwise). 






Fig 3. Winding number problems solved by ecf -modification. 

An example where plain Fourier-based de-Poissonization fails. Rate profiles (a,c,e) 
and imaginary parts of the log characteristic functions {b,d,f) of the "true" model are 
represented by black lines, the empirical analogs by blue or red lines, respectively, if 
the winding number of the initial ecf is zero or non-zero. The deviating rate profiles 
(a) associated with a non-zero winding number of the ecf (fe) are effectively corrected 
by "shrinking" {c,d) and "editing" (e,/), with one exception in case of shrinking. 



A simple remedy is to shrink the loop towards 1 so that it does not circle 
around anymore: given some S e (0, 1), replace 7/i(0) by Jh{0; S) = 6 + {1 — 
d)jhid), and obtain rate estimates Pfc((5) based on 7ft,(-; d) (rather than on 7h(-))- 
As long as 6 is small, the resulting perturbation of the estimates should be small, 
too. For a less crude corrective one starts out directly from the origin of non-zero 
winding numbers: such result if 7/1 considered as a polynomial in w G C has 
zeros within the unit disc (cf. Section 3]). The idea here is to move any such zero 
"outside" by multiplying it with a suitable factor. For definiteness, let e > 
and suppose that jh{w) — X^s^oP'^^'' ~ ^'<i-^i Zi) has zeros ai, . . . ,aK, 
lh{w) = nfcLi('"^ ~ Qf*:)/(1 ~ Q^fe)- Let "edited zeros" 5^ = afc(e) be defined by 
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cife = (1 + e)ak/\ak\ if l^fcl < 1 + and = otherwis^, and put 



k=l 1 - Qffc 

Then define rate estimates Vn{() as above, now replacing 7h(-) by ^h{-\ e). 

Figure [3] illustrates both procedures for the following choice of parameters: 
(j/i, ...,1/4) = (17, 11, 14, 6) and j/^ = for fc > 4; T = 60; /i = 0.050 shrinking 
parameter 5 — 0.02; editing parameter e = 0.075. Our simulations indicate that 
editing induces less bias in the rate estimates than shrinking. Moreover, it works 
always, while for shrinking this holds true only if 5 is selected adaptively. On the 
other hand, such an adaptive choice is easily implemented, and shrinking does 
not require the computation of all zeros of 7^. We will come back to these matters 
elsewhere, along with more extensive simulation studies and applications to 
measured data. 



7. Discussion 

In this paper we have proposed a compound Poisson model along with associated 
statistical procedures for the analysis of higher-order interactions in neuronal 
assembly activity. The proposed model is purely phenomenological; explanation 
of neural mechanisms is neither intended nor within its scope. Its main purpose 
is to provide a formal framework for the statistical analysis of active assemblies 
in sampled data, to the degree they are reflected by the occurrence of higher- 
order spike patterns. This goal is achieved without recourse to simultaneous 
recordings of the single-neuron spike trains. Only the compound signal is used, 
which comprises the superimposed spike trains of all neurons. It is then impossi- 
ble to identify specific groups of neurons that participate in specific coordinated 
events. In return, our approach allows to design statistical procedures that oper- 
ate reliably on much smaller samples than other methods suggested previously. 
Moreover, it avoids certain technical problems associated with the identification 
of individual neurons in population signals typically obtained with extracellu- 
lar recording electrodes: Although spikes have still to be isolated from noise, 
and simultaneous spikes on a single channel have to be correctly identified as 
such, our approach practically eliminates the need of "spike sorting" , a difficult 
and error-prone procedure (often based on a classification of spike wave-forms) 
through which spikes are assigned to individual neurons. 

Technically, our approach hinges on the compound Poisson character of the 
bin count distribution. Ultimately, it is this assumption that establishes the 
connection between the observed counts and the not directly observable higher- 
order correlations in the activity of the neuron population. Doing without any 
such functional relation is impossible in inferences based on the compound ac- 
tivity only. On the other hand, as with any model the assumptions may not be 

^It may be useful to "edit" zeros outside but close to the unit disc as well. Other "zero 
editing" such as reflection at the unit circle may also be worth considering. 

®The critical parameter hvj^ is 2.4 in this case, whereas it was about 1 in the two well- 
behaved examples discussed above. 
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adequate in applications, and potentially serious consequences have to be taken 
into account if they are violated. Let us adress some issues that might turn out 
problematic with actual neuronal behavior. 

It is currently hard to judge how well a real neuronal population conforms to 
the Poisson model, although irregular (i.e. Poissonian) spike trains are consid- 
ered to be characteristic of normal operation of neocortical networks. Recently, 
evidence is accumulating that the single-neuron firing rates in the cortex of be- 
having animals are extremely low (<^; 1 Hz, Lee et al. [31), much lower than 
previously thought (Shadlen and Newsome 26]). In our model, individual spikes 
are distributed over a large number of different assemblies, so "statistics of rare 
events" may be an adequate approximation. Deviations from Poisson statis- 
tics in real neurons, on the other hand, seem to occur mainly for higher firing 
rates, e.g. during stimulus response (Amarashingham et al. 3|) or in the large 
projection neurons of the motor cortex (Baker and Lemon [4|). 

Our assumption of identically distributed bin counts means that the com- 
pound process Z, and the carrier process 1+ , respectively, should have stationary 
increments within the observation period [0, T]. While quasi-stationary stretches 
of data of sufficient length could possibly be obtained in anesthetized animals, 
the task-dependent modulations of neuronal firing rates in behaving animals 
generate inevitable non-stationarities. The consequences of unattended firing 
rate modulations are potentially quite serious: fluctuating rates increase the 
weight in the tail of the count distribution, which falsely predicts the presence 
of (higher-order) correlations. The related problems are lessened in experiments 
avoiding rapid changes of external conditions and long observation periods. In 
regard to the conflict between stationarity and efficiency, which demand short 
and long observation intervals, respectively, the parsimonious parametrization 
of the higher-order interactions appears as a particularly advantageous feature 
of our approach, since it helps to limit the length of the observation period 
required for reliable inferences. 

Alternative parametrizations could also contribute to alleviate the conse- 
quences of instationarity. Thus one could treat the jump rate of the carrier 
process as a nuisance parameter and consider the normalized rates LUn = i'n/i^+ 
as the parameters of interest. Such would be appropriate particularly if only 
the general activation level i'^ is changing across time while the proportions be- 
tween the rates Vn remain unaffected. Stabilizing the variances of estimates may 
also be desirable, which would suggest the parametrization ujn = {vn/v+Y^^ ■ 

For simplicity we assumed that the "coordinated" spikes that constitute a 
pattern in our model of population activity occur strictly simultaneously. Such 
precision of synchrony cannot be expected in reality. Intuitively, this would 
not hamper our theory to the degree to which the "jitter" is small enough 
compared to the width h of the sampling bins, which was 5-20 milliseconds in 
most test simulations we looked at. In fact, this generates yet another conflict of 
interest between "small bins" which improve the performance of the estimators, 
and "large bins" which are less affected by jittered patterns (and, perhaps, by 
residual correlation across time resulting from spike after-effects). 

What could be gained from analyses such as those proposed in this paper? 
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Positive statistical evidence for higher-order correlations in neuronal spike trains 
would be a novel result, potentially challenging conventional rate coding the- 
ories, as well as statements about the sufficiency of second-order correlations 
in capturing "essential aspects" of population activity in certain neuronal sys- 
tems (e.g. Schneidman et al. 25 1, Shlens et al. [13] )■ However, spiking activity 
in a recurrent network will by necessity have a non-trivial statistical structure, 
potentially including also higher-order correlations. This is a consequence of net- 
work features, as for instance shared input and multiple feedback loops. Care- 
ful analysis of corresponding network models is necessary in order to correctly 
identify such structurally determined correlations, and not to mistake them for 
signatures of dynamic neuronal assemblies. 



Appendix A: Validity of the asymptotics 

Our goal here is to investigate under what conditions the heuristic approxima- 
tions of Section [3] can be expected to be valid. Of particular interest are con- 
sistency of estimates and the validity of the first-order (linear) approximation 
(|3.2p . We first discuss a directly tractable special case. 

A.l. A simple test case. Let us consider the asymptotics of the estimate 9+ = 
—h~^ logpo (cf- (|4-^P ) ill some detail. The random variable Lpo is binomially 
distributed with parameters (L,po)- Hence E{po/po) = 1, and Var(po/po) = 
(1 —po)/{Lpq) — > if Lpo — > oo. It follows by a stochastic Taylor expansion of 
the logarithm that under this condition 

r \ 1/2 / r \ 1/2 ^ 
Lpo \ N I Lpo \ ' Po 
h{iy+-u+) = -[- log- 



1 - Po / \^-PoJ Po 



1 -Po/ \Po \ Lpo 

where B — {pq{1—pq)/L) (po — po) is asymptotically standard normally 
distributed if and only if Lpo{l—po) — > oo. We conclude that is asymptotically 
normally distributed, with mean vj^ and variance h^^{\ — pq) / (Lpq) , if and only 
if Lpo{l -po) oo. 

Let us yet explore the consequences of three natural assumptions: 

(a) the intensity of the carrier ppp is bounded away from zero, = 0(1); 

(b) the bin width h stays bounded, h = 0(1); 

(c) P+ is consistent. 

Consistency of P+ means that its asymptotic variance tends to zero, h~^{l — 
Pq)/(Lpo) = (e'"'+ - l)/(T/i) ^ 0. Since v+/T < (0''"^+ - l)/(T/i), this im- 
plies y+jT = 0(1), which is thus a necessary condition for the consistency of 
9+. Incidentally, the /i-independence of the ultimate lower bound v+/T for the 
asymptotic variance shows that choosing the bin width h small is only par- 
tially helpful; given u+ the rate of convergence can be effectively improved only 
by increasing the observation length T. Another implication of consistency is 
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{Lpo{l — po)) ^ — 0(1), i.e., consistency ofD^ implies its asymptotic normality. 
In fact, let us first suppose that hvj^ stays bounded away from zero. Then 



0(1) then by (a) and 
o(l). 



and this tends to zero by (b) and (c). Similarly, if hv^ 
(c) 

hT'h'"'+/{l - e~'"'+) ~ hl{Thv+) = v-"^ {v+lT) 

These considerations, along with the exponential dependence of the asymptotic 
variance on hv^ — — logpoi point to the crucial role of the quantity hv^ — which 
is evident, on the other hand, considering that hvj^ is the expected number of 
points of the carrier process in a bin. The present findings are corroborated in 
the following. 

A. 2. General case. Consistency of the estimates Vn requires that the log- 
ecf, log 7/1, is close in some sense to the logarithm of the common characteristic 
function of the bin counts Zi, log^h- In view of the expansions (j3.ip . (|3.2p . 
closeness largely depends on the order of magnitude of the random variables 
^h{9), which is immediate from (|3.3p : 



= (Th)-' (exp 



2h^^i^k (1 - cos{9k)) 



The last expression is < (Th) ' (e'*'"^+ — 1) for every 9, so on average (across 1 
we have by Jensen's inequality 



> 7^ / E\h-'U0)\'d9 



Th 



2n 

> {Th)-' 



(A.l) 



exp 



^ f 2hY,^''k{l~cos{9k))d9 



- 1 



Th 



Let us remark here that the expression Si := {2Tr)~' J^^ E \h~'£^h{d)\^ d9 rep- 
resents the quantity relevant to uniform consistency of all estimators of the 
form (12. 5|) . This is (i) because to first order, one has for any linear functional 
Ac = J2k '^k the approximation 



Ac Ac ~ ~ — 
27r 



h-'US)C{9)d9- 



and (ii) because the last expression tends to in mean square, uniformly across 
aU Ac such that \\C\\l = (27r)-i/^^ \C{e)\'^ dO < 1, if and only if '~\''^ 0. 
With g{x) = (e^ — ^)/x, the lower bound in (jA.ip becomes 

Si > g{2hv+) 2v+/T > 2v+/T. 
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Consequently, /T — s- is a minimal requirement for the uniform consistency 
of the estimates; it is also sufficient if stays bounded. Summarizing so far, 
the following conditions are found to be necessar'^ for uniform consistency of 
the estimates: 

(CO) T^oo, h^O{l) (so that L = T//i -> oo), and i^+/r 0. 

Going beyond consistency, let us now examine the (uniform) validity of the 
asymptotic (co-)variance approximations. These depend on the validity of the 
linear stochastic approximation (|3.2p , hence on the negligibility of the quadratic 
term (the first one omitted) in p.ip . Its expected mean square — which is 
relevant to the present question for the same reason as Si was for uniform 
consistency — can be estimated as in the following lemma proven below. 

Lemma. Let E2 ^ {2tt)-'^ J^^ E \h-'^ ^,,{0)^1^ d0. Then 

S2<3(^^^^j +3A(e8'^^+-l)^:,,. (A.2) 

Assuming (CO), under what additional conditions is £2 of smaller order than 
the minimum order of the main term. Si? First, suppose that hv+ stays bounded. 
Then g (4/11/4.) ^-nd g{%hvj^^) stay bounded as well, and the main term is of the 
order v^/T. Writing e2 in the form 

^2 = 3/^^ 5(4/..+)) ' + 3 (^) ' ^ g{^h.^) 

one finds that indeed £2 — o{iy+/T) under the assumed conditions. Thus the 
quadratic term is of smaller order than the main term, hence negligible. 
Next, suppose that hi^^ becomes large. Then the main term Si is at least of the 
order e^^''+/{Th), by ((XT|) . while 

,<3(£^) ..r^ + 3^(2^) '.'T^^3(^) (HV^H'T). 

Hence, if 6^'"^+ / (Th) ~ o(l) as is necessary for uniform consistency in this case, 
then the condition h^T = 0(1) implies £2=0 (e^'"'+ /{Th)^ , and the quadratic 
term is negligible compared to the main term. 

Putting everything together, we find that the following conditions additionally 
to (CO) are decisive for the validity of the proposed asymptotics. 

(CI) Either hv+ = 0(1); or hv+ 00, e^'"'+/{Th) 0, and h^T = 0(1). 

Proof of lemma. In order to prove (jA.2p , let us introduce the random variables 
Ci{S) = e*^^'/7(6'). (Henceforth we supress h whenever it appears as a parameter 
only.) Then 

'^Boundedness of h may not really be necessary, but is a natural condition in our setting. 
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and so because of O(^) — 0(" 



E\m'\' = L-'Y: 



E{aAo)-mkM-miA-o)-mi.A-o)-i)- 



Since each single factor in the fourfold product has expectation zero and the Q (9) 
are independent, only those terms contribute to the fourfold sum for which every 
two or all four indices are identical. Therefore, 

= -1 im If ia-9) -if + ^E m ly e (a-e) if 

where ({9) denotes a generic random variable 

Let us now evaluate the single expectations, in reverse order. Similarly as in 
(lO) we get 

E{C{9) - l){C{-9) - 1) =EC{9)a^9) - 1 (A.4) 
= |7(6')r2_i = exp[2/tE'^fe(l-cos(fc6l))] - 1 < 6^'*''+ -1. 

Next, since 

EX EX ^ \EXf < {E\X\f and E {({9) - l){C{-9) - 1) ^ E \(:{9) - if , 



E im - If E {c{-9) If < {E m - If)' 

(exp[2/iE'^fe(l-cos(/c6'))] - l)^ < {e'^'"'+ - if 



(A.5) 



Last, 



E(ae)-if{c{~e)~if 



i0Z 



1 



- 1 



- 1 



1 



hidW 7(0) 7(- 

2 



(A.6) 



1 



1 



1 



1 



17(^)1 



7(2g) 7(-2g) 2 
7(0)2 ^ 7(-0)2 |7(0)P 
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But 



5ft 



7(20) 



< 



7(20) 



exp [ /i ^ i^k (cos(2A:6') - 1) - 2/i ^ t/fe (cos(fc6') - 1) ] 
exp [h^iyk {cos{2ke) - 2 cos(fc6l) + 1) ] 



< phi'+/2 



the latter since cos(2fc6') - 2 cos{k9) + 1 = 2 cos(fc6') (1 - cos{k9)) < 1/2, so that 
in view of (jA.6p we get 



£;(C(0)-i)^(C(-^?)-i)' 

< exp [ 4/i ^ i/fc (1 - cos(fc6i)) ] - 1 + 2{e'"'+^^ - 1) 

< 3(eS'"'+ - 1) . 

Putting everything together we finaUy obtain the bound 



(A.7) 



<3(^^^)\3A(eS-._l) 



and the desired estimate (IA.2[) foUows. 

In passing we note the related lower estimate 



1 /p2/ii/+ _ 1 

— / E \h-^^{0f\'' de > 2(1 - h/T) 



T 



(A. 



It can be established similarly as earlier on noting first that by (jA.3p and (jA.4[) 

i?IC(^)'P > ^%^[i?(C(^?)-i)(C(-e)-i)]' 



L3 
2(i--l) 
L3 



(^exp [ 2/i ^ z^fc (1 - cos(fc0)) ] - 1 



then using the convexity of the function f{x) — (e^ — 1)^ (a; > 0) along with 
Jensen's inequality. 



Appendix B: Evaluation of asymptotic (co-)variances 

All asymptotic (co-)variances of interest to us are of the form J{ni, 712, ai, 0L2) 
for some Uj > and aj € {0, 1} (j — 1, 2), where 



J{ni,n2,ai,a2) 



r/,(0i,02) 



d0i d02 



(27r)2 J_^J_„ (e*^i - l)"i {e^^^ - 1)"^ e*^i"i e*''^'^ 



(B.l) 
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This follows from (|3.6p and the subsequent remark on noting that for an indi- 
vidual rate estimate we have C{9) = e^*^", while for the tail estimate pm we 
have C{0) = e-'^'"/(l - e"*^) = e-^d{m-i) ^ (^^^9 _ particular, 

ascov (9„j , P„ J = T"V(ni,n2,0,0) =T"^fi„j,„2 , (B.2) 
ascov , p,„ J = T"V(rni-l,TO2-l, 1, 1) = I]™i,,„2 , (B.3) 
ascov(pm, ;?„) = T^^ J(m — 1, n, 1, 0) . (B-4) 

To evaluate the double integral (|B.1[) we utilize the explicit representation (|5.ip 
of the kernel r/i(6'i, 02 ), from which it follows that the integrand of (jB.ip depends 
on the arguments 6*^ only via the variables e*^^ {j = 1,2). The substitutions 
6j 1-^ e*^^ = 2;j then give 



If f I / \ dzi dz2 



J(ni,n2,ai,a2) = 77— 2 / / V'ai.aa (^l: ^2) + i (B.5) 



where 

^--(-i'-2) = „ 1).. - 1).. • (B-6) 

No matter what aj G {0,1}, the function ipai,a2{'^iT '^2) is a holomorphic (in 
fact, entire) function of the two complex variables zi, Z2 - Therefore, the bivariate 
version of Cauchy's formula applies, and we obtain the following evaluation of 
(jB.ip (and hence, of the asymptotic covariances). 



. 1 9"^+"^^a„a,(zi,Z2) 

J(ni,n2,ai,a2) = — j r a n-, 



(B.7) 

Zl=Z2=0 



Let us give a few examples. In the test case of Section the asymptotic 
variance of P+ — pi times T was found to be (e'"'+ — l)/h. The same expression 
should obtain for J(0, 0,1,1) = ?/'i,i(0,0), and in fact, it does. For Tasvar(i/i) 
one finds J(l, 1, 0, 0) = e^'^+ (lyi+hi/f), and similar formulae result for Vn , n > 2. 
Examples of covariances are 

ill. 2 = rascov(i?i, 1)2) — e'"^+ hvi (1^2 — 1^1 hi^i 



2 



Tascov(pi, i?i) = e + j^i 



Higher-order formulae are a mess, but can be obtained using a computer algebra 
system. 
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Added Note 

The following two articles relevant to the contents of the present paper came to 
the attention of the authors immediately before publication: 

Buchmann, B. and Griibel, R. (2003). Decompounding: an estimation prob- 
lem for Poisson random sums. Ann. Statist. 31, 1054-1074. 

van Es, B., Gugushvili, S. and Spreij, P. (2007). A kernel type nonparametric 
density estimator for decompounding. Bernoulli 13, 672-694. 
See also the references cited in these works. 
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